########################################################################################################
##                                                                                                    ##
##  This is the replication file for the simulations in Tsai, Tsung-han. 2014.                        ##
##  "A Bayesian Approach to Dynamic Panel Data Models with Rarely Changing Variables."                ##
##                                                                                                    ##
########################################################################################################
####################################################################################
##  CHECK THE BIVARIATE DISTRIBUTIONS IN THE EXPERIMENT WHERE UNIT EFFECTS ARE DRAWN FROM A UNIFORM DISTRIBUTION
cor.zu <- 0.8;

uniform.unieff <- unif.var(n.state=n.state, t=n.time, mu.delta=0, sd.delta=1, cor.xu=0.3, cor.zu=cor.zu);

uniform.unieff$xu.corr; uniform.unieff$zu.corr;

mu.x3 <- uniform.unieff$mu.x3; mu.z3 <- uniform.unieff$mu.z3;
delta <- uniform.unieff$delta;

biv.xu <- kde2d(x= delta, y= mu.x3, n=50);

## PLOT THE BIVARIATE DISTRIBUTION OF x3 (z3) AND UNIT EFFECTS
persp(biv.xu, phi=30, theta=45, shade=0.1, border=NA, main="Bivariate Distribution", xlab="Unit Effects", ylab="Country-Level X3", zlab="Density");

#
plot(x= delta, y= mu.x3, main="Bivariate Distribution", xlab="Unit Effects", ylab="Country-Level X3", type="p", pch="+", col="darkgreen");
contour(biv.xu, add=TRUE);

########################################################################################################
##  Robustness check for the bivariate distribution.                                                  ##
##  The correlation between unit effects and W is varied.                                             ##
##  The correlation between unit effects and X is fixed to 0.3; J=20.                                 ##
##  T=15.                                                                                             ##
########################################################################################################
## Before running, change the working directory in line 50 to the appropriate one.
##
## Require the following files:
## (1) panel_functions.R: auxiliary functions
## (2) bayes.mlm.R: The syntax for the Bayesian multilevel model
## (3) bayes.endo2.R: The syntax for the Bayesian simultaneous equation model (BSEM)
##
## Require the following packages:
## (1) ecodist
## (2) plm
## (3) arm
## (4) R2jags

## CLEAN UP
rm(list=ls());

## SET WORKING DIRECTORY
setwd("");

## LOAD REQUIRED FUNCTIONS
source("panel_functions.R");





## LOAD PACKAGES.
library(ecodist);  # GENERATE CORRELATED SAMPLES: corgen
library(plm);
library(arm);
library(R2jags);

## SET THE SEED, WHICH TELLS US WHERE THE RANDOM PROCESS STARTS.
set.seed(04102014);

## THE NUMBER OF UNITS AND TIME PERIODS
n.state <- 20; n.time <- 15;

####################################################################################
##  RUN SIMULATIONS WHERE UNIT EFFECTS ARE UNIFORMLY DISTRIBUTED IN THE INTERVAL BETWEEN -2 AND 2

## TRUE VALUES OF PARAMETERS
beta <- c(0.8, 0.5, 2, -1.5, -2.5, 1.8, 3);

## GENERATE EXOGENOUS VARIABLES
exog.var <- uncor.var(n.state=n.state, t=n.time, mu.x=c(1,0), mu.z=c(2,2), sd.x=c(3,1), sd.z=c(1,2) );
x1 <- exog.var$x1; x2 <- exog.var$x2; z1 <- exog.var$z1; z2 <- exog.var$z2;

## FINAL STORAGES
corr.out <- list()

ols.est <- list();
fe.est <- list();
re.est <- list();
fevd.est <- list();
mlr.est <- list();
mlm.est <- list();
bmlm.est <- list();
bayes.est <- list();


## THE NUMBER OF REPLICATIONS
n.iter <- 100;


## SET UP THE CORRELATION BETWEEN THE UNIT EFFECTS AND REGRESSORS: x3 OR z3 (THE OTHER IS FIXED AT 0.3)
#cor.xu <- 0.8;

r <- seq(0, 0.9, 0.1);  # EXPERIMENTS FOR DIFFERENT CORRELATION PARAMETERS

for (q in 1:length(r)) {  # LOOP FOR DIFFERENT EXPERIMENTS
cor.zu <- r[q];


## STORAGE FOR ESTIMATES OF CORRELATIONS
CORR <- matrix(NA, ncol=6, nrow=n.iter);
colnames(CORR) <- c("tcor.dx","tcor.dz","rho.dx","rho.dz","sigma.d","sigma.y");

## STORAGES FOR COEFFICIENTS ESTIMATES
coef.ols.matrix <- matrix(NA, nrow=n.iter, ncol=7);
colnames(coef.ols.matrix) <- c("ly","x1","x2","x3","z1","z2","z3");

coef.fe.matrix <- matrix(NA, nrow=n.iter, ncol=7);
colnames(coef.fe.matrix) <- c("ly","x1","x2","x3","z1","z2","z3");

coef.re.matrix <- matrix(NA, nrow=n.iter, ncol=7);
colnames(coef.re.matrix) <- c("ly","x1","x2","x3","z1","z2","z3");

coef.fevd.matrix <- matrix(NA, nrow=n.iter, ncol=7);
colnames(coef.fevd.matrix) <- c("ly","x1","x2","x3","z1","z2","z3");

coef.mlr.matrix <- matrix(NA, nrow=n.iter, ncol=7);
colnames(coef.mlr.matrix) <- c("ly","x1","x2","x3","z1","z2","z3");

coef.mlm.matrix <- matrix(NA, nrow=n.iter, ncol=7);
colnames(coef.mlm.matrix) <- c("ly","x1","x2","x3","z1","z2","z3");

coef.bmlm.matrix <- matrix(NA, nrow=n.iter, ncol=7);
colnames(coef.bmlm.matrix) <- c("ly","x1","x2","x3","z1","z2","z3");

coef.bayes.matrix <- matrix(NA, nrow=n.iter, ncol=7);
colnames(coef.bayes.matrix) <- c("ly","x1","x2","x3","z1","z2","z3");


for (g in 1:n.iter) {  # LOOP FOR REPLICATIONS


## GENERATE CORRELATED VARIABLES
endog.var <- cor.var(n.state=n.state, t=n.time, mu.delta=0, sd.delta=1, cor.xu=0.3, cor.zu= cor.zu);

x3 <- endog.var$x3; z3 <- endog.var$z3
xu.corr <- endog.var$xu.corr; zu.corr <- endog.var$zu.corr;
delta <- endog.var$delta;

## GENERATE THE OUTCOME VARIABLE AND ARRANGE THE DATASET
out.data <- sim.data(n.state=n.state, t=n.time, mu=1, phi=0.8, beta=c(0.5,2,-1.5,-2.5,1.8,3), x1=x1, x2=x2, x3=x3, z1=z1, z2=z2, z3=z3, delta=delta);



########################################
## SIMULATED DATA ANALYSIS: 
########################################
## METHOD 1: POOLED OLS
pooled.ols <- lm(y ~ ly+x1+x2+x3+z1+z2+z3, data=out.data);
#summary(pooled.ols)$coef[,1:2];

coef.ols <- summary(pooled.ols)$coef[-1,1];

## MATRIX OF ESTIMATES
coef.ols.matrix[g,] <- coef.ols;


########################################
## METHOD 2: FIXED EFFECTS MODEL
fe.mod <- plm(y ~ ly+x1+x2+x3+z1+z2+z3, data=out.data, index=c("country", "year"), model="within"); 
#summary(fe.mod)$coef[,1:2];

coef.fe <- summary(fe.mod)$coef[,1];

## MATRIX OF ESTIMATES
coef.fe.matrix[g,] <- coef.fe;


########################################
## METHOD 3: RANDOM EFFECTS MODEL
re.mod <- plm(y ~ ly+x1+x2+x3+z1+z2+z3, data=out.data, index=c("country", "year"), model="random", random.method="amemiya"); 
#summary(re.mod)$coef[,1:2];

coef.re <- summary(re.mod)$coef[-1,1];

## MATRIX OF ESTIMATES
coef.re.matrix[g,] <- coef.re;


########################################
## METHOD 4: FIXED EFFECTS VECTOR DECOMPOSITION (fevd) MODEL
# FIRST STAGE: FIT A STANDARD FIXED EFFECTS MODEL
fevd.1 <- plm(y ~ ly+x1+x2+x3+z1+z2+z3, data=out.data, index=c("country", "year"), model="within"); 
#summary(fevd.1)$coef[,1:2];

b <- coef(fevd.1); b <- matrix(b, nrow=length(b), ncol=1, byrow=FALSE);
X.bar <- out.data[, c(12:18)];
X.bar <- as.matrix(X.bar);
e <- residuals(fevd.1);

e <- matrix(e, nrow=n.state, ncol=n.time, byrow=TRUE);
ebar <- rep(rowMeans(e), each=n.time);
ebar <- matrix(ebar, nrow=n.state*n.time, ncol=1, byrow=FALSE)

u.hat <- out.data$ybar - X.bar %*% b - ebar;

# SECOND STAGE: REGRESS THE ESTIMATED UNIT EFFECTS ON THE SLOWLY CHANGING VARIABLES
fevd.2 <- lm(u.hat ~ -1+out.data$z1+out.data$z2+out.data$z3);
#summary(fevd.2)$coef[,c(1,2)];

h <- residuals(fevd.2);

# THIRD STAGE: RERUN THE FULL MODEL WITHOUT UNIT EFFECTS BUT INCLUDE h
fevd.3 <- lm(out.data$y ~ out.data$ly + out.data$x1 + out.data$x2 + out.data$x3 + out.data$z1 + out.data$z2 + out.data$z3 + h);
#summary(fevd.3)$coef[2:8,1];

coef.fevd <- summary(fevd.3)$coef[2:8,1];

## MATRIX OF ESTIMATES
coef.fevd.matrix[g,] <- coef.fevd;


########################################
## METHOD 5: MULTILEVEL MODEL BY MLE WITHOUT WITHIN-GROUP MEANS
mlm.ind <- lmer(y ~ ly+x1+x2+x3+z1+z2+z3+(1|country), data=out.data);
#summary(mlm.ind);

coef.mlm <- coef(mlm.ind)$country[1,2:8];

## MATRIX OF ESTIMATES
coef.mlr.matrix[g,] <- as.numeric(coef.mlm);


########################################
## METHOD 6: MULTILEVEL MODELING WITH GROUP MEANS
mlm.mean <- lmer(y ~ ly+x1+x2+x3+z1+z2+z3+x3bar+z3bar+(1|country), data=out.data);
#summary(mlm.ind);

coef.mlmean <- coef(mlm.mean)$country[1,2:8];

## MATRIX OF ESTIMATES
coef.mlm.matrix[g,] <- as.numeric(coef.mlmean);


########################################
## METHOD 7: BAYESIAN MULTILEVEL MODEL WITH GROUP MEANS
N <- dim(out.data)[1];
J <- length(unique(out.data$country));
country <- out.data$country

y <- out.data$y

# COVARIATES OF Y REGRESSION
X.1 <- cbind(out.data$ly, out.data$x1, out.data$x2, out.data$x3, out.data$z1, out.data$z2, out.data$z3);
M.1 <- dim(X.1)[2];

# COVARIATES OF THE SECOND LEVEL REGRESSION
Z.1 <- cbind(1, unique(out.data$x3bar), unique(out.data$z3bar) );
L <- dim(Z.1)[2];

data1 <- list("y", "N", "J", "country", "X.1", "M.1", "Z.1", "L");

# INITIAL VALUES
inits1 <- list(
	list( tau.y=1, B=rep(0, M.1), tau.d=0.1, G=rep(-3, L) ),
	list( tau.y=0.1, B=rep(3, M.1), tau.d=0.5, G=rep(3, L) ),
	list( tau.y=0.5, B=rep(-3, M.1), tau.d=1, G=rep(0, L) )
	)

parameters1 <- c("B");

n.chains <- length(inits1)
n.iteration <- 10000;
frac.discard <- 0.5;  # FRACTION OF SAMPLES DISCARDED
n.discard <- n.iteration*frac.discard; n.discard;  # BURNIN-IN PERIOD
n.thining <- 5;
n.saved.chin <- (n.iteration - n.discard)/n.thining; n.saved.chin;  # SAMPLES AFTER THINING SAVED FOR EACH CHAIN
n.saved.chin*n.chains;  # TOTAL SAMPLES SAVED


## hand-off to JAGS




bmlm <- jags(model.file="bayes.mlm.R",
             data = data1,
             inits = inits1,
             parameters.to.save = parameters1,
             n.chains = n.chains,   # (DEFAULT: 1)
             n.iter = n.iteration,
             n.burnin = n.discard,
             n.thin = n.thining,
             DIC = TRUE
             );

#print(bmlm);

## 
mcmcburnin.mcmclist <- as.mcmc(bmlm);

	
# POSTERIOR STATS
mcmc.burnin <- as.matrix(mcmcburnin.mcmclist);
std <- apply(mcmc.burnin, 2, sd);
Mean <- apply(mcmc.burnin, 2, mean);
coef.bmlm <- cbind(Mean[1:7], std[1:7]);

## MATRIX OF ESTIMATES
coef.bmlm.matrix[g,] <- coef.bmlm[,1];


########################################
## METHOD 8: BAYESIAN SIMULTANEOUS EQUATION MODEL

# re-run lines 244-248

# COVARIATES OF Y REGRESSION
X.1 <- cbind(1, out.data$ly, out.data$x1, out.data$x2, out.data$x3, out.data$z1, out.data$z2, out.data$z3);
M.1 <- dim(X.1)[2];

# COVARIATES OF X3 REGRESSION
X.2 <- matrix(rep(1,N),ncol=1)
L <- dim(X.2)[2];

# COVARIATES OF W3 REGRESSION
W.1 <- matrix(rep(1,N),ncol=1)
Q <- dim(W.1)[2];

K <- 3 # THE NUMBER OF REGRESSIONS

W <- diag(K)

data1 <- list("y", "N", "J", "country", "X.1", "M.1", "X.2", "L", "W.1", "Q", "K", "W");


# INITIAL VALUES
inits1 <- list(
	list(tau.y=1, tau.x3=1, tau.z3=1, B=rep(0, M.1), G=rep(0, L), H=rep(0, Q), D=array(rep(0,J*K), c(J,K)), Tau.D=W ),
	list(tau.y=.1, tau.x3=.1, tau.z3=.1, B=rep(-3, M.1), G=rep(-3, L), H=rep(-3, Q), D=array(rep(3,J*K), c(J,K)), Tau.D=W ),
	list(tau.y=.5, tau.x3=.5, tau.z3=.5, B=rep(3, M.1), G=rep(3, L), H=rep(3, Q), D=array(rep(-3,J*K), c(J,K)), Tau.D=W )
	)

parameters1 <- c("sigma.y", "B", "sigma.d", "sigma.x", "rho.dx", "rho.dz", "sigma.z")

n.chains <- length(inits1)
n.iteration <- 10000;
frac.discard <- 0.5;  # FRACTION OF SAMPLES DISCARDED
n.discard <- n.iteration*frac.discard; n.discard;  # BURNIN-IN PERIOD
n.thining <- 5;
n.saved.chin <- (n.iteration - n.discard)/n.thining; n.saved.chin  # SAMPLES AFTER THINING SAVED FOR EACH CHAIN
n.saved.chin*n.chains  # TOTAL SAMPLES SAVED


## hand-off to JAGS




bml1 <- jags(model.file="bayes.endo2.R",
             data = data1,
             inits = inits1,
             parameters.to.save = parameters1,
             n.chains = n.chains,   # (DEFAULT: 1)
             n.iter = n.iteration,
             n.burnin = n.discard,
             n.thin = n.thining,
             DIC = TRUE
             );

#print(bml1);


## 
mcmcburnin.mcmclist <- as.mcmc(bml1);








	
# POSTERIOR STATS
mcmc.burnin <- as.matrix(mcmcburnin.mcmclist);
#dim(mcmc.burnin);
#Mode <- apply(mcmc.burnin, 2, quantile, prob=c(0.025,0.5,0.975));
std <- apply(mcmc.burnin, 2, sd);
Mean <- apply(mcmc.burnin, 2, mean);
coef.bayes <- cbind(Mean[2:8], std[2:8]);

## MATRIX OF ESTIMATES
coef.bayes.matrix[g,] <- coef.bayes[,1];


A <- c()
A[1] <- endog.var$xu.corr; 
A[2] <- endog.var$zu.corr; 
A[3:6] <- Mean[c("rho.dx","rho.dz","sigma.d","sigma.y")]

CORR[g,] <- A;


}  # END OF REPLICATIONS



########################################
##  SAVE THE RESULTS INTO THE FINAL STORAGES
corr.out[[q]] <- CORR;

ols.est[[q]] <- coef.ols.matrix;
fe.est[[q]] <- coef.fe.matrix;
re.est[[q]] <- coef.re.matrix;
fevd.est[[q]] <- coef.fevd.matrix;
mlr.est[[q]] <- coef.mlr.matrix;
mlm.est[[q]] <- coef.mlm.matrix;
bmlm.est[[q]] <- coef.bmlm.matrix;
bayes.est[[q]] <- coef.bayes.matrix;


} # END OF EACH EXPERIMENT



##
length(ols.est);
dim(ols.est[[1]]);

ols.est;
fe.est;
re.est;
fevd.est;
mlr.est;
mlm.est;
bmlm.est
bayes.est;


##  SAVE THE OUTPUT OF EXPERIMENTS
setwd("/robustness2/");
write.csv(corr.out, "sim_ro2_corr.csv", row.names=FALSE);

write.csv(ols.est, "est_ro2_ols.csv", row.names=FALSE);
write.csv(fe.est, "est_ro2_fe.csv", row.names=FALSE);
write.csv(re.est, "est_ro2_re.csv", row.names=FALSE);
write.csv(fevd.est, "est_ro2_fevd.csv", row.names=FALSE);
write.csv(mlr.est, "est_ro2_mlr.csv", row.names=FALSE);
write.csv(mlm.est, "est_ro2_mlm.csv", row.names=FALSE);
write.csv(bmlm.est, "est_ro2_bmlm.csv", row.names=FALSE);
write.csv(bayes.est, "est_ro2_bayes.csv", row.names=FALSE);




################################################################################
# REPORTING RESULTS OF SIMULATIONS
################################################################################
###################################################################################
##  Experiment 2: The correlation between unit effects and W is varied.                               
##  The correlation between unit effects and X is fixed to 0.3.                                      
###################################################################################

## CLEAN UP
rm(list=ls());

##
est.mc <- read.csv("/robustness2/est_ro2_ols.csv");

## TRUE VALUES OF PARAMETERS
beta <- c(0.8, 0.5, 2, -1.5, -2.5, 1.8, 3);

beta.v <- rep(beta, times=10);
beta.m <- matrix(rep(beta.v, times=100), nrow=100, ncol=70, byrow=TRUE);

##  CALCULATE AVERAGED ABSOLUTE BIAS OVER TEN EXPERIMENTS
bias.est <- abs(est.mc - beta.m);

bias.m <- matrix(rep(0, 100*7), nrow=100, ncol=7);

for (i in 1:10) {
	bias.m <- bias.m + bias.est[,(1+7*(i-1)):(7+7*(i-1))];
}

colnames(bias.m) <- c("ly","x1","x2","x3","z1","z2","z3");

bias.m <- bias.m/10;

##  CALCULATE THE MEAN OF AVERAGED ABSOLUTE BIAS ACROSS 500 REPLICATIONS
bias.avg <- apply(bias.m, 2, mean);

##  CALCULATE MONTE CARLO ERRORS (MCE) OF AVERAGED ABSOLUTE BIAS
bias.mce <- apply(bias.m, 2, sd);


##  CALCULATE RMSE OVER TEN EXPERIMENTS
bias.sq <- (est.mc - beta.m)^2;

bias.sq.m <- matrix(rep(0, 100*7), nrow=100, ncol=7);

for (i in 1:10) {
	bias.sq.m <- bias.sq.m + bias.sq[,(1+7*(i-1)):(7+7*(i-1))];
}

colnames(bias.sq.m) <- c("ly","x1","x2","x3","z1","z2","z3");

bias.sq.m <- sqrt(bias.sq.m/10);


##  CALCULATE THE MEAN OF RMSE ACROSS 500 REPLICATIONS
rmse.avg <- apply(bias.sq.m, 2, mean);

##  CALCULATE MONTE CARLO ERRORS (MCE) OF RMSE
rmse.mce <- apply(bias.sq.m, 2, sd);


##
round(rbind(bias.avg, bias.mce), 3)[,c("ly", "x3", "z3")];

round(rbind(rmse.avg, rmse.mce), 3)[,c("ly", "x3", "z3")];


# OLS AAVB
            ly    x3    z3
bias.avg 0.019 0.138 0.176
bias.mce 0.003 0.021 0.028

# OLS RMSE
rmse.avg 0.022 0.157 0.208
rmse.mce 0.003 0.022 0.030


##  FOR FIGURE 1
ols.bias <- matrix(apply(bias.est, 2, mean), nrow=10, ncol=7, byrow=TRUE);
colnames(ols.bias) <- c("ly","x1","x2","x3","z1","z2","z3");
rownames(ols.bias) <- seq(0, 0.9, 0.1);

ols.rmse <- matrix(sqrt(apply(bias.sq, 2, mean)), nrow=10, ncol=7, byrow=TRUE);
colnames(ols.rmse) <- c("ly","x1","x2","x3","z1","z2","z3");
rownames(ols.rmse) <- seq(0, 0.9, 0.1);


##
est.mc <- read.csv("/robustness2/est_ro2_fe.csv");

##  RE-RUN LINES 475-518

# FE AAVB
            ly    x3   z3
bias.avg 0.006 0.050 0.122
bias.mce 0.001 0.011 0.032

# FE RMSE
rmse.avg 0.007 0.061 0.150
rmse.mce 0.002 0.014 0.038


##  FOR FIGURE 
fe.bias <- matrix(apply(bias.est, 2, mean), nrow=10, ncol=7, byrow=TRUE);
colnames(fe.bias) <- c("ly","x1","x2","x3","z1","z2","z3");
rownames(fe.bias) <- seq(0, 0.9, 0.1);

fe.rmse <- matrix(sqrt(apply(bias.sq, 2, mean)), nrow=10, ncol=7, byrow=TRUE);
colnames(fe.rmse) <- c("ly","x1","x2","x3","z1","z2","z3");
rownames(fe.rmse) <- seq(0, 0.9, 0.1);


##
est.mc <- read.csv("/robustness2/est_ro2_fevd.csv");

##  RE-RUN LINES 475-518

# FEVD AAVB
            ly    x3    z3
bias.avg 0.006 0.050 0.359
bias.mce 0.001 0.011 0.046

# FEVD RMSE
rmse.avg 0.007 0.061 0.403
rmse.mce 0.002 0.014 0.049


##  FOR FIGURE 
fevd.bias <- matrix(apply(bias.est, 2, mean), nrow=10, ncol=7, byrow=TRUE);
colnames(fevd.bias) <- c("ly","x1","x2","x3","z1","z2","z3");
rownames(fevd.bias) <- seq(0, 0.9, 0.1);

fevd.rmse <- matrix(sqrt(apply(bias.sq, 2, mean)), nrow=10, ncol=7, byrow=TRUE);
colnames(fevd.rmse) <- c("ly","x1","x2","x3","z1","z2","z3");
rownames(fevd.rmse) <- seq(0, 0.9, 0.1);


##
est.mc <- read.csv("/robustness2/est_ro2_re.csv");

##  RE-RUN LINES 475-518

# RE AAVB
            ly    x3    z3
bias.avg 0.005 0.048 0.128
bias.mce 0.001 0.011 0.027

# RE RMSE
rmse.avg 0.006 0.060 0.155
rmse.mce 0.002 0.013 0.030


##  FOR FIGURE 
re.bias <- matrix(apply(bias.est, 2, mean), nrow=10, ncol=7, byrow=TRUE);
colnames(re.bias) <- c("ly","x1","x2","x3","z1","z2","z3");
rownames(re.bias) <- seq(0, 0.9, 0.1);

re.rmse <- matrix(sqrt(apply(bias.sq, 2, mean)), nrow=10, ncol=7, byrow=TRUE);
colnames(re.rmse) <- c("ly","x1","x2","x3","z1","z2","z3");
rownames(re.rmse) <- seq(0, 0.9, 0.1);


##
est.mc <- read.csv("/robustness2/est_ro2_mlm.csv");

##  RE-RUN LINES 475-518

# MLM WITH GROUP MEANS, AAVB
            ly    x3    z3
bias.avg 0.005 0.050 0.118
bias.mce 0.001 0.011 0.032

# MLM WITH GROUP MEANS, RMSE
rmse.avg 0.006 0.061 0.145
rmse.mce 0.002 0.013 0.037


##  FOR FIGURE 
mlm.bias <- matrix(apply(bias.est, 2, mean), nrow=10, ncol=7, byrow=TRUE);
colnames(mlm.bias) <- c("ly","x1","x2","x3","z1","z2","z3");
rownames(mlm.bias) <- seq(0, 0.9, 0.1);

mlm.rmse <- matrix(sqrt(apply(bias.sq, 2, mean)), nrow=10, ncol=7, byrow=TRUE);
colnames(mlm.rmse) <- c("ly","x1","x2","x3","z1","z2","z3");
rownames(mlm.rmse) <- seq(0, 0.9, 0.1);


##
est.mc <- read.csv("/robustness2/est_ro2_bayes.csv");

##  RE-RUN LINES 475-518

# BAYESIAN AAVB
            ly    x3    z3
bias.avg 0.005 0.049 0.113
bias.mce 0.001 0.011 0.029

# BAYESIAN RMSE
rmse.avg 0.006 0.060 0.140
rmse.mce 0.002 0.013 0.033


##  FOR FIGURE 
bayes.bias <- matrix(apply(bias.est, 2, mean), nrow=10, ncol=7, byrow=TRUE);
colnames(bayes.bias) <- c("ly","x1","x2","x3","z1","z2","z3");
rownames(bayes.bias) <- seq(0, 0.9, 0.1);

bayes.rmse <- matrix(sqrt(apply(bias.sq, 2, mean)), nrow=10, ncol=7, byrow=TRUE);
colnames(bayes.rmse) <- c("ly","x1","x2","x3","z1","z2","z3");
rownames(bayes.rmse) <- seq(0, 0.9, 0.1);



###################  MAKE PLOTS
par(mfrow=c(3,2),mar=c(0,0,2,0),oma=c(4,4,2,4));

######################## ABSOLUTE BIAS OF ly FROM xcor
corr <- seq(0, 0.9, 0.1);

plot(ols.bias[,1] ~ corr, type="n", xlab="", ylab="Absolute Bias", main="(a) Absolute Bias of LDV", xlim=c(0,1.1), ylim=c(0, 0.03), yaxt="n", xaxt="n");

axis(side=2, las=1, tick=TRUE, cex.axis=1);

lines(corr, ols.bias[,1], lty=1, col="black", lwd=2);
lines(corr, fe.bias[,1], lty=2, col="red", lwd=2);
lines(corr, fevd.bias[,1], lty=3, col="forestgreen", lwd=2);
lines(corr, re.bias[,1], lty=4, col="goldenrod4", lwd=2);
lines(corr, mlm.bias[,1], lty=5, col="darkslategray", lwd=2);
lines(corr, bayes.bias[,1], lty=6, col="darkblue", lwd=2);

text.loc <- cbind(c(0.3, 0.97, 0.97, 0.96), c(0.02, 0.007, 0.0055, 0.004));
ml.name <- c("pooled OLS", "FE, FEVD", "RE", "MLM, BSEM");
text(text.loc, ml.name, cex=0.8, col="black");

#mtext("Absolute Bias", side=2, line=2.5, outer=F, cex=.8);


######################## RMSE OF ly FROM xcor
plot(ols.rmse[,1] ~ corr, type="n", xlab="", ylab="RMSE", main="(b) RMSE of LDV", xlim=c(0,1.1), ylim=c(0, 0.03), yaxt="n", xaxt="n");

axis(side=4, las=1, tick=TRUE, cex.axis=1);

lines(corr, ols.rmse[,1], lty=1, col="black", lwd=2);
lines(corr, fe.rmse[,1], lty=2, col="red", lwd=2);
lines(corr, fevd.rmse[,1], lty=3, col="forestgreen", lwd=2);
lines(corr, re.rmse[,1], lty=4, col="goldenrod4", lwd=2);
lines(corr, mlm.rmse[,1], lty=5, col="darkslategray", lwd=2);
lines(corr, bayes.rmse[,1], lty=6, col="darkblue", lwd=2);

text.loc <- cbind(c(0.3, 0.97, 0.97, 0.96), c(0.023, 0.009, 0.007, 0.005));
ml.name <- c("pooled OLS", "FE, FEVD", "RE", "MLM, BSEM");
text(text.loc, ml.name, cex=0.8, col="black");

#mtext("RMSE", side=4, line=2.5, outer=F, cex=1);


######################## ABSOLUTE BIAS OF x3 FROM xcor
plot(ols.bias[,4] ~ corr, type="n", xlab="", ylab="Absolute Bias", main="(c) Absolute Bias of x3", xlim=c(0,1.1), ylim=c(0, 0.45), yaxt="n", xaxt="n");

axis(side=2, las=1, tick=TRUE, cex.axis=1);

lines(corr, ols.bias[,4], lty=1, col="black", lwd=2);
lines(corr, fe.bias[,4], lty=2, col="red", lwd=2);
lines(corr, fevd.bias[,4], lty=3, col="forestgreen", lwd=2);
lines(corr, re.bias[,4], lty=4, col="goldenrod4", lwd=2);
lines(corr, mlm.bias[,4], lty=5, col="darkslategray", lwd=2);
lines(corr, bayes.bias[,4], lty=6, col="darkblue", lwd=2);

text.loc <- cbind(c(0.9, 0.97, 0.9), c(0.43, 0.2, 0.05));
ml.name <- c("pooled OLS", "RE", "FE, FEVD, MLM, BSEM");
text(text.loc, ml.name, cex=0.8, col="black");

#mtext("Absolute Bias", side=2, line=2.5, outer=F, cex=.8)


######################## RMSE OF x3 FROM xcor
plot(ols.rmse[,4] ~ corr, type="n", xlab="", ylab="RMSE", main="(d) RMSE of x3", xlim=c(0,1.1), ylim=c(0, 0.45), yaxt="n", xaxt="n");

axis(side=4, las=1, tick=TRUE, cex.axis=1);

lines(corr, ols.rmse[,4], lty=1, col="black", lwd=2);
lines(corr, fe.rmse[,4], lty=2, col="red", lwd=2);
lines(corr, fevd.rmse[,4], lty=3, col="forestgreen", lwd=2);
lines(corr, re.rmse[,4], lty=4, col="goldenrod4", lwd=2);
lines(corr, mlm.rmse[,4], lty=5, col="darkslategray", lwd=2);
lines(corr, bayes.rmse[,4], lty=6, col="darkblue", lwd=2);

text.loc <- cbind(c(0.92, 0.97, 0.9), c(0.43, 0.2, 0.05));
ml.name <- c("pooled OLS", "RE", "FE, FEVD, MLM, BSEM");
text(text.loc, ml.name, cex=0.8, col="black");

#mtext("RMSE", side=4, line=2.5, outer=F, cex=1)


######################## ABSOLUTE BIAS OF w3 FROM xcor
plot(ols.bias[,7] ~ corr, type="n", xlab="", ylab="Absolute Bias", main="(e) Absolute Bias of w3", xlim=c(0,1.1), ylim=c(0, 0.4), yaxt="n", xaxt="n");

axis(side=2, las=1, tick=TRUE, cex.axis=1);

lines(corr, ols.bias[,7], lty=1, col="black", lwd=2);
lines(corr, fe.bias[,7], lty=2, col="red", lwd=2);
lines(corr, fevd.bias[,7], lty=3, col="forestgreen", lwd=2);
lines(corr, re.bias[,7], lty=4, col="goldenrod4", lwd=2);
lines(corr, mlm.bias[,7], lty=5, col="darkslategray", lwd=2);
lines(corr, bayes.bias[,7], lty=6, col="darkblue", lwd=2);

text.loc <- cbind(c(0.97, 0.99, 0.97, 0.99), c(0.07, 0.12, 0.32, 0.1));
ml.name <- c("pooled OLS", "FE, MLM", "FEVD", "BSEM, RE");
text(text.loc, ml.name, cex=0.8, col="black");

axis(side=1, las=1, tick=TRUE, cex.axis=1);

#mtext("Absolute Bias", side=2, line=2.5, outer=F, cex=.8)


######################## RMSE OF w3 FROM xcor
plot(ols.rmse[,7] ~ corr, type="n", xlab="", ylab="RMSE", main="(f) RMSE of w3", xlim=c(0,1.1), ylim=c(0, 0.4), yaxt="n", xaxt="n");

axis(side=4, las=1, tick=TRUE, cex.axis=1);

lines(corr, ols.rmse[,7], lty=1, col="black", lwd=2);
lines(corr, fe.rmse[,7], lty=2, col="red", lwd=2);
lines(corr, fevd.rmse[,7], lty=3, col="forestgreen", lwd=2);
lines(corr, re.rmse[,7], lty=4, col="goldenrod4", lwd=2);
lines(corr, mlm.rmse[,7], lty=5, col="darkslategray", lwd=2);
lines(corr, bayes.rmse[,7], lty=6, col="darkblue", lwd=2);

text.loc <- cbind(c(0.97, 0.99, 0.97, 0.99), c(0.1, 0.15, 0.36, 0.13));
ml.name <- c("pooled OLS", "FE, MLM", "FEVD", "BSEM, RE");
text(text.loc, ml.name, cex=0.8, col="black");

axis(side=1, las=1, tick=TRUE, cex.axis=1);

#mtext("RMSE", side=4, line=2.5, outer=F, cex=1)

mtext(expression(paste("corr(x3,",delta,")")), side=1, line=2.5, outer=T, cex=1);


###################################################################################
## RESULTS OF CORRELATION
###################################################################################
## CLEAN UP
rm(list=ls());

## 
xcor_corr <- read.csv("/robustness2/sim_ro2_corr.csv");
dim(xcor_corr);

xcor_corr[1,1:6];

avg.corr <- colMeans(xcor_corr);
mavg.corr <- matrix(avg.corr, ncol=6, byrow=TRUE);
colnames(mavg.corr) <- c("tcor.dx","tcor.dz","rho.dx","rho.dz","sigma.d","sigma.y");
rownames(mavg.corr) <- seq(0,0.9,0.1);
round(mavg.corr, 3);

    tcor.dx tcor.dz rho.dx rho.dz sigma.d sigma.y
0     0.220   0.006  0.192 -0.007   1.014   1.000
0.1   0.214   0.085  0.216  0.024   1.015   1.010
0.2   0.219   0.175  0.179  0.158   0.967   1.001
0.3   0.220   0.270  0.203  0.176   1.011   0.999
0.4   0.214   0.359  0.196  0.298   1.003   1.002
0.5   0.222   0.438  0.198  0.362   1.013   1.017
0.6   0.226   0.539  0.216  0.452   0.966   1.009
0.7   0.217   0.651  0.196  0.523   1.030   1.004
0.8   0.217   0.751  0.191  0.624   1.006   1.010
0.9   0.220   0.868  0.168  0.682   1.001   1.001


# SINCE THE DISTRIBUTIONS ARE NOT SYMMETRIC, HIGHEST POSTERIOR DENSITY (HPD) INTERVALS ARE CALCULATED
## USE "HPDint" FUNCTION
## LOAD REQUIRED FUNCTIONS
source("panel_functions.R");
hpdint <- HPDint(xcor_corr, prob=0.95);

hpdint <- t(hpdint$HPDinterval);

corr.med <- apply(xcor_corr, 2, quantile, probs=0.5);

corr.hpd <- rbind(hpdint, corr.med);

corr.com <- rbind(corr.hpd[1,], corr.hpd[3,], corr.hpd[2,]);

#corr.com <- apply(xcor_corr, 2, quantile, probs=c(0.1,0.5,0.9));
sub <- seq(1,60,by=6);

t.cor.dx <- t(corr.com[,seq(1,60,by=6)]);
t.cor.dz <- t(corr.com[,seq(2,60,by=6)]);
cor.dx <- t(corr.com[,seq(3,60,by=6)]);
cor.dz <- t(corr.com[,seq(4,60,by=6)]);

r <- seq(0,0.9,0.1);


##  FIGURE fr21
plot(y=t.cor.dx[,2], x=r, ylim=c(-.4, .9), type="p", pch="+", ylab="", xlab=expression(paste("corr(x3,",delta,")")), col="black", lwd=4, main="", axes=FALSE, cex.lab=1.5);

segments(0, seq(-0.9,0.9,0.1), 1, seq(-0.9,0.9,0.1), col="grey70", lty=2, lwd=0.5);

segments(r, t.cor.dx[,1], r, t.cor.dx[,3], col="grey30", lwd=3);
points(y=t.cor.dx[,2], x=r, pch="+", cex=1);
axis(side=1, at=seq(0, .9, .1), tick=TRUE, cex.axis=1.2);
axis(side=2, at=seq(-1, 1, .1), tick=TRUE, cex.axis=1.2, las=1);


#lines(x=r, y=cor.dx[,2], type="b", pch="+", col="blue", lwd=4);
segments(r+0.01, cor.dx[,1], r+0.01, cor.dx[,3], col="black", lwd=3, lty=2);
points(y=cor.dx[,2], x=r+0.01, pch="+", cex=1.5, col="black");


##  FIGURE fr22
plot(y=t.cor.dz[,2], x=r, ylim=c(-.4, .9), type="p", pch="+", ylab="", xlab=expression(paste("corr(w3,",delta,")")), col="black", lwd=4, main="", axes=FALSE, cex.lab=1.5);

segments(0, seq(-0.9,0.9,0.1), 1, seq(-0.9,0.9,0.1), col="grey70", lty=2, lwd=0.5);

segments(r, t.cor.dz[,1], r, t.cor.dz[,3], col="grey30", lwd=3);
points(y=t.cor.dz[,2], x=r, pch="+", cex=1);
axis(side=1, at=seq(0, .9, .1), tick=TRUE, cex.axis=1.2);
axis(side=2, at=seq(-1, 1, .1), tick=TRUE, cex.axis=1.2, las=1);


#lines(x=r, y=cor.dz[,2], type="b", pch="+", col="blue", lwd=4);
segments(r+0.01, cor.dz[,1], r+0.01, cor.dz[,3], col="black", lwd=3, lty=2);
points(y=cor.dz[,2], x=r+0.01, pch="+", cex=1.5, col="black");




